
clear
rng default


%% Empirical distribution from data on cohort 1940-1941
%%Men%%
year=1940;
% Probability distribution of Ychosen|X
condprobsmen_temp=load('probs_allyears_data.mat',['condprobsmen' num2str(year) '_140']);
condprobsmen=condprobsmen_temp.(['condprobsmen' num2str(year) '_140']);
PYX_cond=[condprobsmen(:,2:end) condprobsmen(:,1)]; %5x6
%Order: P(Ychosen=1|X=1) P(Ychosen=2|X=1) P(Ychosen=3|X=1) P(Ychosen=4|X=1) P(Ychosen=5|X=1) P(Ychosen=0|X=1) 
%       P(Ychosen=1|X=2) P(Ychosen=2|X=2) P(Ychosen=3|X=2) P(Ychosen=4|X=2) P(Ychosen=5|X=2) P(Ychosen=0|X=2) 
%       P(Ychosen=1|X=3) P(Ychosen=2|X=3) P(Ychosen=3|X=3) P(Ychosen=4|X=3) P(Ychosen=5|X=3) P(Ychosen=0|X=3) 
%       P(Ychosen=1|X=4) P(Ychosen=2|X=4) P(Ychosen=3|X=4) P(Ychosen=4|X=4) P(Ychosen=5|X=4) P(Ychosen=0|X=4) 
%       P(Ychosen=1|X=5) P(Ychosen=2|X=5) P(Ychosen=3|X=5) P(Ychosen=4|X=5) P(Ychosen=5|X=5) P(Ychosen=0|X=5) 
probsmen_temp=load('probs_allyears_data.mat',['probsmen' num2str(year) '_140']);
PX=probsmen_temp.(['probsmen' num2str(year) '_140']); %5x1
%Order: P(X=1) P(X=2) P(X=3) P(X=4) P(X=5) 

%Regroup to have 2 types plus singles
PYX_joint=PYX_cond.*kron(PX, ones(1,6));
PYX_joint_new_1=[PYX_joint(:,1)+PYX_joint(:,2)  PYX_joint(:,3)+PYX_joint(:,4)+PYX_joint(:,5)  PYX_joint(:,6)];
PYX_joint_new=[PYX_joint_new_1(1,:)+PYX_joint_new_1(2,:); PYX_joint_new_1(3,:)+PYX_joint_new_1(4,:)+PYX_joint_new_1(5,:)]; 
PX_emp=[PX(1)+PX(2); PX(3)+PX(4)+PX(5)];
PYX_cond_emp=PYX_joint_new./(kron(PX_emp, ones(1,3)));
PYX_cond_emp=[PYX_cond_emp(1,:) PYX_cond_emp(2,:)];


%%Women%%
year=1941;
% Probability distribution of Xchosen|Y
condprobswomen_temp=load('probs_allyears_data.mat',['condprobswomen' num2str(year) '_140']);
condprobswomen=condprobswomen_temp.(['condprobswomen' num2str(year) '_140']);
PXY_cond=[condprobswomen(:,2:end) condprobswomen(:,1)]; %5x6
%Order: P(Xchosen=1|Y=1) P(Xchosen=2|Y=1) P(Xchosen=3|Y=1) P(Xchosen=4|Y=1) P(Xchosen=5|Y=1) P(Xchosen=0|Y=1) 
%       P(Xchosen=1|Y=2) P(Xchosen=2|Y=2) P(Xchosen=3|Y=2) P(Xchosen=4|Y=2) P(Xchosen=5|Y=2) P(Xchosen=0|Y=2) 
%       P(Xchosen=1|Y=3) P(Xchosen=2|Y=3) P(Xchosen=3|Y=3) P(Xchosen=4|Y=3) P(Xchosen=5|Y=3) P(Xchosen=0|Y=3) 
%       P(Xchosen=1|Y=4) P(Xchosen=2|Y=4) P(Xchosen=3|Y=4) P(Xchosen=4|Y=4) P(Xchosen=5|Y=4) P(Xchosen=0|Y=4) 
%       P(Xchosen=1|Y=5) P(Xchosen=2|Y=5) P(Xchosen=3|Y=5) P(Xchosen=4|Y=5) P(Xchosen=5|Y=5) P(Xchosen=0|Y=5) 
probswomen_temp=load('probs_allyears_data.mat',['probswomen' num2str(year) '_140']);
PY=probswomen_temp.(['probswomen' num2str(year) '_140']); %5x1
%Order: P(Y=1) P(Y=2) P(Y=3) P(Y=4) P(Y=5) 

%Regroup to have 2 types plus singles
PXY_joint=PXY_cond.*kron(PY, ones(1,6));
PXY_joint_new_1=[PXY_joint(:,1)+PXY_joint(:,2)  PXY_joint(:,3)+PXY_joint(:,4)+PXY_joint(:,5)  PXY_joint(:,6)];
PXY_joint_new=[PXY_joint_new_1(1,:)+PXY_joint_new_1(2,:); PXY_joint_new_1(3,:)+PXY_joint_new_1(4,:)+PXY_joint_new_1(5,:)]; 
PY_emp=[PY(1)+PY(2); PY(3)+PY(4)+PY(5)];
PXY_cond_emp=PXY_joint_new./(kron(PY_emp, ones(1,3)));
PXY_cond_emp=[PXY_cond_emp(1,:) PXY_cond_emp(2,:)];


%% Get the systematic match surplus
n_mixture=2;
ntypes_m=2;
ntypes_w=2;

variance_mixture=linspace(1,50,n_mixture); 
corr_mixture=linspace(1,-10,n_mixture);
MU_mixture=zeros(n_mixture, ntypes_w+1);
SIGMA_mixture=zeros(ntypes_w+1,ntypes_w+1,n_mixture);
for j=1:n_mixture
    SIGMA_mixture(:,:,j)=[variance_mixture(j) corr_mixture(j) corr_mixture(j); ...
                          corr_mixture(j) variance_mixture(j) corr_mixture(j); ...
                          corr_mixture(j) corr_mixture(j) variance_mixture(j)];                
end
gm = gmdistribution(MU_mixture, SIGMA_mixture);
gw = gm;


number_starting=10^4; %number of starting values for fminunc
r=10^5; %number of simulations to evaluate expectation of the max
%Simulations to evaluate expectation of the max
epsilon_sim=random(gm,r); %rxntypes_w+1 (epsilon1, epsilon2, ..., epsilon0) 
eta_sim=random(gw,r); %rxntypes_m+1 (eta1, eta2, ..., eta0)
%Options for fimunc
options =  optimoptions(@fminunc, 'MaxFunctionEvaluations', 10^4 ,'MaxIterations', 10^4, 'StepTolerance', 10^(-8), 'Display', 'off');
%Matrix of starting values for fimunc
u_starting=normrnd(0,1,number_starting,ntypes_m); %rxntypes_m (u1x u2x)
v_starting=normrnd(0,1,number_starting,ntypes_w); %rxntypes_m (u1x u2x)

G_star2_u11 = @(x) G_star(x, PYX_cond_emp(2),epsilon_sim, u_starting, number_starting,options); %derivative wrto P(Y=1|X=1)
G_star2_u21 = @(x) G_star(PYX_cond_emp(1),x,epsilon_sim, u_starting, number_starting,options); %derivative wrto P(Y=2|X=1)
G_star2_u12 = @(x) G_star(x, PYX_cond_emp(5),epsilon_sim, u_starting, number_starting,options); %derivative wrto P(Y=1|X=2)
G_star2_u22 = @(x) G_star(PYX_cond_emp(4),x,epsilon_sim, u_starting, number_starting,options); %derivative wrto P(Y=2|X=2)

G_star2_v11 = @(x) G_star(x, PXY_cond_emp(2),eta_sim, v_starting, number_starting,options); %derivative wrto P(X=1|Y=1)
G_star2_v21 = @(x) G_star(PXY_cond_emp(1),x,eta_sim, v_starting, number_starting,options); %derivative wrto P(X=2|Y=1)
G_star2_v12 = @(x) G_star(x, PXY_cond_emp(5),eta_sim, v_starting, number_starting,options); %derivative wrto P(X=1|Y=2)
G_star2_v22 = @(x) G_star(PXY_cond_emp(4),x,eta_sim, v_starting, number_starting,options); %derivative wrto P(X=2|Y=2)

h=10^(-4);
u11_emp=(G_star2_u11(PYX_cond_emp(1)+h)-G_star2_u11(PYX_cond_emp(1)))/h;
u21_emp=(G_star2_u21(PYX_cond_emp(2)+h)-G_star2_u21(PYX_cond_emp(2)))/h;
u12_emp=(G_star2_u12(PYX_cond_emp(4)+h)-G_star2_u12(PYX_cond_emp(4)))/h;
u22_emp=(G_star2_u22(PYX_cond_emp(5)+h)-G_star2_u22(PYX_cond_emp(5)))/h;
v11_emp=(G_star2_v11(PXY_cond_emp(1)+h)-G_star2_v11(PXY_cond_emp(1)))/h;
v21_emp=(G_star2_v21(PXY_cond_emp(2)+h)-G_star2_v21(PXY_cond_emp(2)))/h;
v12_emp=(G_star2_v12(PXY_cond_emp(4)+h)-G_star2_v12(PXY_cond_emp(4)))/h;
v22_emp=(G_star2_v22(PXY_cond_emp(5)+h)-G_star2_v22(PXY_cond_emp(5)))/h;


Phi11_emp=u11_emp+v11_emp;
Phi12_emp=u21_emp+v12_emp;
Phi21_emp=u12_emp+v21_emp;
Phi22_emp=u22_emp+v22_emp;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% Primitives
nm=6000; %number of men
nw=6000; %number of women
supX=[1 2]; %support of X (men's types)
supY=[1 2]; %support of Y (women's types)
supXY=sup(supX, supY); %support of (X,Y)
%supXY order:
%(1,1)
%(2,1)
%(1,2)
%(2,2)
PhiXY=[Phi11_emp Phi12_emp; Phi21_emp Phi22_emp];
%PhiXY=
%Phi11 Phi12
%Phi21 Phi22



%% Generate data
%Men
X=discrete(nm, round(PX_emp,3), supX); %nmx1; 
epsilon=epsilon_sim(1:nm,:);
%epsilon=random(gm,nm);%nmxsize(supX,2)+1;
%epsilon=(epsilon1, epsilon2, ..., epsilon0) from Gaussian mixture, independent of X

%Women
Y=discrete(nw, round(PY_emp,3), supY); %nwx1;     
eta=eta_sim(1:nw,:);
%eta = random(gw,nw); %nwxsize(supY,2)+1; 
%eta=(eta1, eta2, ..., eta0) from Gaussian mixture, independent of Y

%% Equilibrium
%(1) Compute total surplus for each possible couple
Phi=zeros(nm,nw); %nmxnw
for j=1:nw %for each woman
    for i=1:nm %for each man
        Phi(i,j)=PhiXY(X(i), Y(j))+epsilon(i,Y(j))+eta(j,X(i));
    end
end

%(2) Reshape Phi into a vector (nm*nw)x1 following the order
%man1 woman1
%man2 woman1
%...
%man1 woman2
%man2 woman2
%...
%man1 woman3
%man2 woman3
%...
Phir=reshape(Phi,nm*nw,1); %(nm*nw)x1

%(3) Augment Phir for singles following the order
%man1 woman1
%man1 woman2
%man1 woman3
%...
%man2 woman1
%man2 woman2
%man2 woman3
%...
%man1 single
%man2 single
%...
%woman1 single
%woman2 single
%woman3 single
%...

Phira=[Phir; epsilon(:,end); eta(:,end)]; %(nm*nw+nm+nw)x1

%(3) Find equilibrium matching 
mu=equilibrium(nw, nm,Phira); %(nm*nw+nm+nw)x1


%% Conditional probabilities
% Men
Ychosen=zeros(nm,1); %nmx1 listing the type of woman man i is matched with
for i=1:nm
    %select from mu the (nw+1)x1 sub-vector reporting the choice of man i
    index=(1:nm:(nm*nw+1)).'+(i-1);
    submu=mu(index);
    %save the type of the woman man i is matched with
    if submu(end)==0 %man i has not chosen to be single
        Ychosen(i)=Y(logical(submu(1:end-1)));
    else %man i has chosen to be single
        Ychosen(i)=0;
    end
end

% Probability distribution of Ychosen|X
PX=[sum(X==1)/nm sum(X==2)/nm];
PYX=[sum(X==1 & Ychosen==1)/nm sum(X==1 & Ychosen==2)/nm sum(X==1 & Ychosen==0)/nm ...
     sum(X==2 & Ychosen==1)/nm sum(X==2 & Ychosen==2)/nm sum(X==2 & Ychosen==0)/nm];
PYX_cond=PYX./kron(PX,ones(1,size(supY,2)+1)); 
%order: P(Ychosen=1|X=1) P(Ychosen=2|X=1) P(Ychosen=0|X=1) P(Ychosen=1|X=2) P(Ychosen=2|X=2) P(Ychosen=0|X=2)


% Women
Xchosen=zeros(nw,1); %nwx1 listing the type of man woman j is matched with
for j=1:nw
    %select from mu the (nm+1)x1 sub-vector reporting the choice of woman j
    index=[(1:nm).'+(j-1)*nm; nm*nw+nm + j];
    submu=mu(index);
    %save the type of the man woman j is matched with
    if submu(end)==0 %woman j has not chosen to be single
        Xchosen(j)=X(logical(submu(1:end-1)));
    else %woman j has chosen to be single
        Xchosen(j)=0;
    end
end

% Probability distribution of Xchosen|Y
PY=[sum(Y==1)/nw sum(Y==2)/nw];
PXY=[sum(Y==1 & Xchosen==1)/nw sum(Y==1 & Xchosen==2)/nw sum(Y==1 & Xchosen==0)/nw ...
     sum(Y==2 & Xchosen==1)/nw sum(Y==2 & Xchosen==2)/nw sum(Y==2 & Xchosen==0)/nw];
PXY_cond=PXY./kron(PY,ones(1,size(supX,2)+1)); 
%order: P(Xchosen=1|Y=1) P(Xchosen=2|Y=1) P(Xchosen=0|Y=1) P(Xchosen=1|Y=2) P(Xchosen=2|Y=2) P(Xchosen=0|Y=2)

save('PYX_cond.mat', 'PYX_cond')
save('PXY_cond.mat', 'PXY_cond')

%% Systematic payoff - Logit assumption
%Notation: u12=man is type 2 and chooses woman of type 1
%u11=u1(1), u12=u1(2), u22=u2(2), u21=u2(1);
u11=log(PYX_cond(1)/PYX_cond(3)); %u1(1)=log(P(Ychosen=1|X=1)/P(Ychosen=0|X=1))
u12=log(PYX_cond(4)/PYX_cond(6)); %u1(2)=log(P(Ychosen=1|X=2)/P(Ychosen=0|X=2))
u22=log(PYX_cond(5)/PYX_cond(6)); %u2(2)=log(P(Ychosen=2|X=2)/P(Ychosen=0|X=2))
u21=log(PYX_cond(2)/PYX_cond(3)); %u2(1)=log(P(Ychosen=2|X=1)/P(Ychosen=0|X=1))
%Notation: v12=woman is type 2 and chooses man of type 1
%v11=v1(1), v12=v1(2), v22=v2(2), v21=v2(1); 
v11=log(PXY_cond(1)/PXY_cond(3)); %v1(1)=log(P(Xchosen=1|Y=1)/P(Xchosen=0|Y=1))
v12=log(PXY_cond(4)/PXY_cond(6)); %v1(2)=log(P(Xchosen=1|Y=2)/P(Xchosen=0|Y=2))
v22=log(PXY_cond(5)/PXY_cond(6)); %v2(2)=log(P(Xchosen=2|Y=2)/P(Xchosen=0|Y=2))
v21=log(PXY_cond(2)/PXY_cond(3)); %v2(1)=log(P(Xchosen=2|Y=1)/P(Xchosen=0|Y=1))          


U_11_CS=u11;
U_12_CS=u21;
U_21_CS=u12;
U_22_CS=u22;
save('U_CS.mat', 'U_11_CS', 'U_12_CS', 'U_21_CS', 'U_22_CS')

V_11_CS=v11;
V_12_CS=v12;
V_21_CS=v21;
V_22_CS=v22;
save('V_CS.mat', 'V_11_CS', 'V_12_CS', 'V_21_CS', 'V_22_CS')

%% Systematic payoff - Truth 

G_star2_u11 = @(x) G_star(x, PYX_cond(2),epsilon_sim, u_starting, number_starting,options); %derivative wrto P(Y=1|X=1)
G_star2_u21 = @(x) G_star(PYX_cond(1),x,epsilon_sim, u_starting, number_starting,options); %derivative wrto P(Y=2|X=1)
G_star2_u12 = @(x) G_star(x, PYX_cond(5),epsilon_sim, u_starting, number_starting,options); %derivative wrto P(Y=1|X=2)
G_star2_u22 = @(x) G_star(PYX_cond(4),x,epsilon_sim, u_starting, number_starting,options); %derivative wrto P(Y=2|X=2)

G_star2_v11 = @(x) G_star(x, PXY_cond(2),eta_sim, v_starting, number_starting,options); %derivative wrto P(X=1|Y=1)
G_star2_v21 = @(x) G_star(PXY_cond(1),x,eta_sim, v_starting, number_starting,options); %derivative wrto P(X=2|Y=1)
G_star2_v12 = @(x) G_star(x, PXY_cond(5),eta_sim, v_starting, number_starting,options); %derivative wrto P(X=1|Y=2)
G_star2_v22 = @(x) G_star(PXY_cond(4),x,eta_sim, v_starting, number_starting,options); %derivative wrto P(X=2|Y=2)


u11_true=(G_star2_u11(PYX_cond(1)+h)-G_star2_u11(PYX_cond(1)))/h;
u21_true=(G_star2_u21(PYX_cond(2)+h)-G_star2_u21(PYX_cond(2)))/h;
u12_true=(G_star2_u12(PYX_cond(4)+h)-G_star2_u12(PYX_cond(4)))/h;
u22_true=(G_star2_u22(PYX_cond(5)+h)-G_star2_u22(PYX_cond(5)))/h;
v11_true=(G_star2_v11(PXY_cond(1)+h)-G_star2_v11(PXY_cond(1)))/h;
v21_true=(G_star2_v21(PXY_cond(2)+h)-G_star2_v21(PXY_cond(2)))/h;
v12_true=(G_star2_v12(PXY_cond(4)+h)-G_star2_v12(PXY_cond(4)))/h;
v22_true=(G_star2_v22(PXY_cond(5)+h)-G_star2_v22(PXY_cond(5)))/h;


%% Scale normalisations - Truth
%If we do not impose independence
U_11_true=u11_true/(u11_true/u11); %=U_11_CS
U_12_true=u21_true/(u11_true/u11); %different from U_12_CS
U_21_true=u12_true/(u12_true/u12);%=U_21_CS
U_22_true=u22_true/(u12_true/u12);%different from U_22_CS
save('U_true.mat', 'U_11_true', 'U_12_true', 'U_21_true', 'U_22_true')

V_11_true=v11_true/(v11_true/v11); %=V_11_CS
V_21_true=v21_true/(v11_true/v11); %different from V_21_CS
V_12_true=v12_true/(v12_true/v12); %=V_12_CS
V_22_true=v22_true/(v12_true/v12); %different from V_22_CS
save('V_true.mat', 'V_11_true', 'V_12_true', 'V_21_true', 'V_22_true')

%If we impose independence
U_11_true_indep=u11_true/(u11_true/u11); %=U_11_CS=U_11_true
U_12_true_indep=u21_true/(u11_true/u11); %=U_12_true
U_21_true_indep=u12_true/(u11_true/u11);
U_22_true_indep=u22_true/(u11_true/u11);
save('U_true_indep.mat', 'U_11_true_indep', 'U_12_true_indep', 'U_21_true_indep', 'U_22_true_indep')

V_11_true_indep=v11_true/(v11_true/v11); %=V_11_CS=V_11_true
V_21_true_indep=v21_true/(v11_true/v11); %=V_21_true
V_12_true_indep=v12_true/(v11_true/v11);
V_22_true_indep=v22_true/(v11_true/v11);
save('V_true_indep.mat', 'V_11_true_indep', 'V_12_true_indep', 'V_21_true_indep', 'V_22_true_indep')


